!!****m* ABINIT/m_melemts
!! NAME
!!  m_melemts
!!
!! FUNCTION
!!  This module defines an object used as database to store matrix
!!  elements of several potentials and operators between two Bloch states.
!!  These values are used in the GW part of abinit to evaluate QP energies
!!  using the perturbative approach.
!!
!! COPYRIGHT
!! Copyright (C) 2008-2020 ABINIT group (MG)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!!
!! NOTES
!!  * This module is supposed to be used only in the GW part to facilitate
!!    further developments. The object might change in the future. Thus
!!    contact Matteo Giantomassi if you wish to use this piece of code
!!    for your developments. In particular we might decide to switch
!!    to ragged arrays Mels(nkibz,nsppol*nspinor**2)%data
!!
!!  * Routines tagged with "@type_name" are tightly connected to the definition of the data type.
!!    Tightly connected means that the proper functioning of the implementation relies on the
!!    assumption that the tagged procedure is consistent with the type declaration.
!!    Every time a developer changes the structure "type_name" adding new entries, he/she has to make sure
!!    that all the tightly connected routines are changed accordingly to accommodate the modification of the data type.
!!    Typical examples of tightly connected routines are creation, destruction or reset methods.
!!
!! TODO
!!  This module can be moved to a higher level directory.
!!
!! PARENTS
!!
!! CHILDREN
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

#include "abi_common.h"

MODULE m_melemts

 use defs_basis
 use m_errors
 use m_xmpi
 use m_abicore

 use m_fstrings,       only : tolower, sjoin
 use m_numeric_tools,  only : print_arr

 implicit none

 private
!!***

!----------------------------------------------------------------------

!!****t* m_melemts/melflags_t
!! NAME
!!
!! FUNCTION
!!  Container for the flags defining the status of the corresponding
!!  pointer defined in the type melements_t. Possible values are:
!!    * 0 if the correspondending array is not allocated.
!!    * 1 if allocated but not yet calculated.
!!    * 2 if allocated and calculated
!!
!! SOURCE

 type,public :: melflags_t

  integer :: has_hbare=0
  integer :: has_lexexch=0
  integer :: has_sxcore=0
  integer :: has_vhartree=0
  integer :: has_vu=0
  integer :: has_vxc=0
  integer :: has_vxcval=0
  integer :: has_vxcval_hybrid=0
  integer :: only_diago=0
  ! 1 if only diagonal elements are calculated

 end type melflags_t

 public :: melflags_reset  ! Reset the value of the flags.
 public :: melflags_copy   ! Copy the object
!!***

!----------------------------------------------------------------------

!!****t* m_melemts/melements_t
!! NAME
!!
!! FUNCTION
!!  Structure defining a database to store the matrix elements of operators
!!  needed for GW calculations.
!!
!! SOURCE

 type,public :: melements_t

  integer :: nkibz
  ! Number of k-points.

  integer :: nsppol
  ! Number of independent spin-polarizations.

  integer :: nspinor
  ! 1 for collinear, 2 for noncollinear.

  integer :: nspden
  ! Number of independent spin-density components.

  integer :: bmin,bmax
  ! min and Max band index over k-points and spin.
  ! Used to dimension the arrays below.

  integer, allocatable :: bands_idx(:,:,:)
  ! bands_idx(2,nkibz,nsppol)
  ! min and Max band index for each k-point and spin.

  integer, allocatable :: iscalc(:,:)
  ! stat_k(nkibz,nsppol)
  ! 1 if this k-point and spin has been calculated, 0 otherwise.

  real(dp), allocatable :: kibz(:,:)
  ! kibz(3,nkibz)
  ! The list of k-points in reduced coordinates.

  complex(dpc), allocatable :: hbare(:,:,:,:)
  ! hbare(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the bare Hamiltonian.

  complex(dpc), allocatable :: sxcore(:,:,:,:)
  ! sxcore(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the Fock operator generated by core electrons.

  complex(dpc), allocatable :: vhartree(:,:,:,:)
  ! vhartree(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the Hartree potential.

  complex(dpc), allocatable :: vlexx(:,:,:,:)
  ! vlexx(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the local exact exchange potential.

  complex(dpc), allocatable :: vu(:,:,:,:)
  ! vu(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the U Hamiltonian.

  complex(dpc), allocatable :: vxc(:,:,:,:)
  ! vxc(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of XC potential, including core if present.

  complex(dpc), allocatable :: vxcval(:,:,:,:)
  ! vxcval(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the XC potential, valence-only contribution.

  complex(dpc), allocatable :: vxcval_hybrid(:,:,:,:)
  ! vxcval_hybrid(b1:b2,b1:b2,nkibz,nsppol*nspinor**2)
  ! Matrix elements of the XC potential for hybrid calculations, valence-only contribution.

  type(melflags_t) :: flags

 end type melements_t

 public :: melements_init      ! Initialize the object
 public :: melements_free      ! Free memory
 public :: melements_herm      ! Construct the lower triangle from the upper triangle
 public :: melements_mpisum    ! Perform a collective SUM within the MPI communicator comm
 public :: melements_print     ! Print matrix elements
 public :: melements_zero      ! Set matrix elements connecting states with different irrep to zero.
 !public :: mels_get_exene_core
!!***

 integer,parameter,private :: NNAMES=8
 integer,parameter,private :: NAMELEN=13

! List of matrix element names.
! Important note: Remember to add a new entry here when you add a new matrix element.
! then change the routine my_select_melements accordingly.
 character(len=NAMELEN),parameter,private :: &
& ANAMES(NNAMES) = (/"vxc          ",&
&                    "vxcval       ",&
&                    "vxcval_hybrid",&
&                    "sxcore       ",&
&                    "vu           ",&
&                    "vlexx        ",&
&                    "vhartree     ",&
&                    "hbare        "/)

CONTAINS  !========================================================================================

!!****f* m_melemts/melflags_reset
!! NAME
!! melflags_reset
!!
!! FUNCTION
!!  Set all flags in melflags_t to 0.
!!
!! INPUTS
!!
!! PARENTS
!!      m_melemts,sigma
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine melflags_reset(Mflags)

!Arguments ------------------------------------
 type(melflags_t),intent(inout) :: Mflags

! *************************************************************************

 ! @melflags_t
 Mflags%has_hbare           = 0
 Mflags%has_sxcore          = 0
 Mflags%has_vhartree        = 0
 Mflags%has_vu              = 0
 Mflags%has_vxc             = 0
 Mflags%has_vxcval          = 0
 Mflags%has_vxcval_hybrid   = 0
 Mflags%has_lexexch         = 0
 Mflags%only_diago          = 0

end subroutine melflags_reset
!!***

!----------------------------------------------------------------------

!!****f* m_melemts/melflags_copy
!! NAME
!! melflags_copy
!!
!! FUNCTION
!!  Copy an object storing the flags.
!!
!! INPUTS
!!  Mflags_in=The flags to be copied.
!!
!! OUTPUT
!!  Mflags_out=The new set of flags.
!!
!! PARENTS
!!      m_melemts
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine melflags_copy(Mflags_in, Mflags_out)

!Arguments ------------------------------------
 type(melflags_t),intent(in)    :: Mflags_in
 type(melflags_t),intent(inout) :: Mflags_out

! *************************************************************************

 call melflags_reset(Mflags_out)

 ! @melflags_t
 Mflags_out%has_hbare           = Mflags_in%has_hbare
 Mflags_out%has_sxcore          = Mflags_in%has_sxcore
 Mflags_out%has_vhartree        = Mflags_in%has_vhartree
 Mflags_out%has_vu              = Mflags_in%has_vu
 Mflags_out%has_vxc             = Mflags_in%has_vxc
 Mflags_out%has_vxcval          = Mflags_in%has_vxcval
 Mflags_out%has_vxcval_hybrid   = Mflags_in%has_vxcval_hybrid
 Mflags_out%has_lexexch         = Mflags_in%has_lexexch
 Mflags_out%only_diago          = Mflags_in%only_diago

end subroutine melflags_copy
!!***

!----------------------------------------------------------------------

!!****f* m_melemts/melements_free
!! NAME
!! melements_free
!!
!! FUNCTION
!!  Free all dynamic memory of the database
!!
!! INPUTS
!!  Mels<melements_t>=The database to be freed
!!
!! OUTPUT
!!  See side effects
!!
!! PARENTS
!!      sigma
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine melements_free(Mels)

!Arguments ------------------------------------
 type(melements_t),intent(inout) :: Mels

! *************************************************************************

 DBG_ENTER("COLL")

 ! @melements_t

!integer arrays
 ABI_SFREE(Mels%bands_idx)
 ABI_SFREE(Mels%iscalc)

!real arrays
 ABI_SFREE(Mels%kibz)

!complex arrays
 ABI_SFREE(Mels%hbare)
 ABI_SFREE(Mels%sxcore)
 ABI_SFREE(Mels%vhartree)
 ABI_SFREE(Mels%vlexx)
 ABI_SFREE(Mels%vu)
 ABI_SFREE(Mels%vxc)
 ABI_SFREE(Mels%vxcval)
 ABI_SFREE(Mels%vxcval_hybrid)

 ! Reset all has_* flags.
 call melflags_reset(Mels%flags)

 DBG_EXIT("COLL")

end subroutine melements_free
!!***

!----------------------------------------------------------------------

!!****f* m_melemts/my_select_melements
!! NAME
!!  my_select_melements
!!
!! FUNCTION
!!  Helper function returning a pointer to the array "aname" as well as the status of the array
!!
!! INPUTS
!!  Mels<melements_t>=The database.
!!  aname=String with the name of the array.
!!
!! OUTPUT
!!  flag_p=Pointer to the integer defining the status of the array, see melflags_t.
!!  arr_p=The pointer to the array.
!!
!! PARENTS
!!      m_melemts
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine my_select_melements(Mels,aname,flag_p,arr_p)

!Arguments ------------------------------------
!scalars
 integer,pointer :: flag_p
 character(len=*),intent(in) :: aname
 type(melements_t),target,intent(in) :: Mels
!arrays
 complex(dpc),ABI_CONTIGUOUS pointer :: arr_p(:,:,:,:)

! *************************************************************************

 SELECT CASE (tolower(aname))
 CASE ("hbare")
   flag_p => Mels%flags%has_hbare
   arr_p  => Mels%hbare
 CASE ("sxcore")
   flag_p => Mels%flags%has_sxcore
   arr_p  => Mels%sxcore
 CASE ("vhartree")
   flag_p => Mels%flags%has_vhartree
   arr_p  => Mels%vhartree
 CASE ("vlexx")
   flag_p => Mels%flags%has_lexexch
   arr_p  => Mels%vlexx
 CASE ("vu")
   flag_p => Mels%flags%has_vu
   arr_p  => Mels%vu
 CASE ("vxc")
   flag_p => Mels%flags%has_vxc
   arr_p  => Mels%vxc
 CASE ("vxcval")
   flag_p => Mels%flags%has_vxcval
   arr_p  => Mels%vxcval
 CASE ("vxcval_hybrid")
   flag_p => Mels%flags%has_vxcval_hybrid
   arr_p  => Mels%vxcval_hybrid
 CASE DEFAULT
   MSG_ERROR(sjoin("Wrong aname: ", aname))
 END SELECT

end subroutine my_select_melements
!!***

!----------------------------------------------------------------------

!!****f* m_melemts/melements_init
!! NAME
!! melements_init
!!
!! FUNCTION
!!  Initialize the database, allocate arrays according to
!!  Mflags_in, zeroing the content of the allocated arrays.
!!
!! INPUTS
!!  nspinor=Number of spinor components
!!  nsppol=Number of independent spin polarizations.
!!  nspden=Number of spin-density components
!!  b1,b2=min and Max band index over spin and k-points.
!!
!! OUTPUT
!!  Mels=The initialized database with dimensions and allocated memory.
!!
!! PARENTS
!!      calc_vhxc_me
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine melements_init(Mels,Mflags_in,nsppol,nspden,nspinor,nkibz,kibz,bands_idx)

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: nspinor,nspden,nsppol,nkibz
 type(melflags_t),intent(in) :: Mflags_in
 type(melements_t),intent(out) :: Mels
!arrays
 integer,intent(in) :: bands_idx(2,nkibz,nsppol)
 real(dp),intent(in) :: kibz(3,nkibz)

!Local variables-------------------------------
 integer :: ikibz,isppol,bmin,bmax,b1,b2

! *************************************************************************

 DBG_ENTER("COLL")

 !@melements_t

 ! * Copy flags.
 call melflags_copy(Mflags_in, Mels%flags)

 ! * Copy dimensions.
 Mels%nkibz   = nkibz
 Mels%nsppol  = nsppol
 Mels%nspinor = nspinor
 Mels%nspden  = nspden

 ABI_MALLOC(Mels%bands_idx,(2,nkibz,nsppol))
 Mels%bands_idx=bands_idx

 ABI_MALLOC(Mels%iscalc,(nkibz,nsppol))
 Mels%iscalc=0

 bmin = HUGE(1); bmax =-HUGE(1)
 do isppol=1,Mels%nsppol
   do ikibz=1,Mels%nkibz
     if (ANY(Mels%bands_idx(:,ikibz,isppol)/=0)) then
       b1 = Mels%bands_idx(1,ikibz,isppol)
       b2 = Mels%bands_idx(2,ikibz,isppol)
       Mels%iscalc(ikibz,isppol)=1
       bmin = MIN(bmin,b1)
       bmax = MAX(bmax,b2)
       ABI_CHECK(b2>=b1 .and. b1>0,"Wrong b1, b2")
     end if
   end do
 end do

 if (bmin==HUGE(1).or.bmax==-HUGE(1)) then
   MSG_BUG("Wrong bands_idx")
 end if

 Mels%bmin = bmin
 Mels%bmax = bmax

 b1 = Mels%bmin
 b2 = Mels%bmax

! real arrays
 ABI_MALLOC(Mels%kibz,(3,nkibz))
 Mels%kibz = kibz

! complex arrays
 if (Mels%flags%has_hbare == 1) then
   ABI_CALLOC(Mels%hbare,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_sxcore == 1) then
   ABI_CALLOC(Mels%sxcore,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_vhartree == 1) then
   ABI_CALLOC(Mels%vhartree,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_lexexch == 1) then
   ABI_CALLOC(Mels%vlexx,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_vu == 1) then
   ABI_CALLOC(Mels%vu,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_vxc == 1) then
   ABI_CALLOC(Mels%vxc,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_vxcval == 1) then
   ABI_CALLOC(Mels%vxcval,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 if (Mels%flags%has_vxcval_hybrid == 1) then
   ABI_CALLOC(Mels%vxcval_hybrid,(b1:b2,b1:b2,nkibz,nsppol*nspinor**2))
 end if

 DBG_EXIT("COLL")

end subroutine melements_init
!!***

!----------------------------------------------------------------------

!!****f* m_melemts/melements_herm
!! NAME
!! melements_herm
!!
!! FUNCTION
!!  Reconstruc the lower triangle of all calculated arrays.
!!  Assuming Hermitian operator. Works both for collinear and  non-collinear cases.
!!
!! INPUTS
!!  Mels=The database
!!  [aname]=The name of the array to be symmetrized, by default
!!    all calculated arrays are filled.
!!
!! SIDE EFFECTS
!!  All arrays whose flag is 2, are filled assuming an Hermitian operator.
!!
!! PARENTS
!!      calc_vhxc_me
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine melements_herm(Mels,aname)

!Arguments ------------------------------------
!scalars
 character(len=*),optional,intent(in) :: aname
 type(melements_t),intent(inout) :: Mels

!Local variables-------------------------------
 integer :: is,ik,ib,jb,iab,iab_tr,iname
 integer,pointer :: flag_p
 character(len=NAMELEN) :: key
!arrays
 integer,parameter :: trsp_idx(2:4) = [2,4,3]
 complex(dpc),ABI_CONTIGUOUS pointer :: arr_p(:,:,:,:)

! *************************************************************************

 ! === Symmetrize matrix elements ===
 ! * In the collinear case, generate the lower triangle by just doing a complex conjugate.
 ! * In the noncollinear case do also a transposition since A_{12}^{ab} = A_{21}^{ba}^*
 !   2-->2, 3-->4, 4-->3
 !
 do iname=1,NNAMES
   key = ANAMES(iname)
   if (PRESENT(aname)) then
     if (key /= aname) CYCLE
   end if

   call my_select_melements(Mels,key,flag_p,arr_p)

   if (flag_p>0) then
     do ik=1,Mels%nkibz
       do is=1,Mels%nsppol

         do jb=Mels%bmin,Mels%bmax
           do ib=Mels%bmin,jb ! Upper triangle

             if (ib/=jb) then
               arr_p(jb,ib,ik,is)=CONJG(arr_p(ib,jb,ik,is))
               if (Mels%nspinor==2) then
                 do iab=2,4
                   iab_tr=trsp_idx(iab)
                   arr_p(jb,ib,ik,iab)=CONJG(arr_p(ib,jb,ik,iab_tr))
                 end do
               end if
             else ! For ib==jb force real-valued
               arr_p(jb,ib,ik,is)=half*(arr_p(jb,ib,ik,is)+CONJG(arr_p(jb,ib,ik,is)))
               if (Mels%nspinor==2) arr_p(jb,ib,ik,2)=half*(arr_p(ib,jb,ik,2)+CONJG(arr_p(ib,jb,ik,2)))
             end if

           end do !ib
         end do !jb

       end do !is
     end do !ik
   end if

 end do !inames

end subroutine melements_herm
!!***

!----------------------------------------------------------------------

!!****f* m_melemts/melements_mpisum
!! NAME
!! melements_mpisum
!!
!! FUNCTION
!!  Perform a collective SUM within the MPI communicator comm
!!  of the matrix elements stored in the database.
!!
!! INPUTS
!!  Mels=The database
!!  [aname]=The name of a particular array to be summed, by default
!!    all allocated arrays are considered.
!!
!! SIDE EFFECTS
!!  All arrays whose flag==1 are summed within the MPI communicator comm.
!!  In output the corresponding flas is set to 2.
!!
!! PARENTS
!!      calc_vhxc_me
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine melements_mpisum(Mels,comm,aname)

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: comm
 character(len=*),optional,intent(in) :: aname
 type(melements_t),intent(inout) :: Mels

!Local variables-------------------------------
 integer :: iname,ierr
 integer,pointer :: flag_p
 character(len=NAMELEN) :: key
 character(len=500) :: msg
!arrays
 complex(dpc),ABI_CONTIGUOUS pointer :: arr_p(:,:,:,:)

! *************************************************************************

 do iname=1,NNAMES
   key = ANAMES(iname)
   if (PRESENT(aname)) then
    if (key /= aname) CYCLE
   end if

   call my_select_melements(Mels,key,flag_p,arr_p)

   if (flag_p==1) then
     call xmpi_sum(arr_p,comm,ierr)
     if (ierr/=0) then
       write(msg,'(a,i4,2a)')" xmpi_sum reported ierr= ",ierr," for key= ",TRIM(key)
       MSG_ERROR(msg)
     end if
     flag_p=2 ! Tag this array as calculated
   end if
 end do

end subroutine melements_mpisum
!!***

!----------------------------------------------------------------------

!!****f* m_melemts/melements_print
!! NAME
!! melements_print
!!
!! FUNCTION
!!  Printout of the content of all calculated array.
!!  Optionally, it is possible to print the content of a single entry of the database.
!!
!! INPUTS
!!  Mels=The database
!!  [unit]=the unit number for output, defaults to std_out
!!  [prtvol]=verbosity level, defaults to 0
!!  [mode_paral]=either "COLL" or "PERS", default to "COLL"
!!  [header]=title for info
!!
!! OUTPUT
!!  Only writing
!!
!! PARENTS
!!      sigma
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine melements_print(Mels,names_list,header,unit,prtvol,mode_paral)

!Arguments ------------------------------------
!scalars
 type(melements_t),intent(in) :: Mels
 integer,optional,intent(in) :: prtvol,unit
 character(len=*),optional,intent(in) :: names_list(:)
 character(len=*),optional,intent(in) :: header
 character(len=4),optional,intent(in) :: mode_paral

!Local variables-------------------------------
 integer :: my_unt,my_prtvol,max_r,max_c,ii
 integer :: isppol,ikibz,iab,ib,b1,b2,my_nkeys,ikey
 integer,pointer :: flag_p
 character(len=4) :: my_mode
 character(len=NAMELEN) :: key
 character(len=500) :: msg,str,fmt
!arrays
 integer,allocatable :: tab(:)
 character(len=NAMELEN),allocatable :: my_keys(:)
 complex(dpc),allocatable :: mat(:,:)

 type rarr_dpc4
   complex(dpc),ABI_CONTIGUOUS pointer :: arr_p(:,:,:,:)
 end type rarr_dpc4
 type(rarr_dpc4),allocatable :: data_p(:)

! *************************************************************************

 !@melements_t
 my_unt   =std_out; if (PRESENT(unit      )) my_unt   =unit
 my_prtvol=0      ; if (PRESENT(prtvol    )) my_prtvol=prtvol
 my_mode  ='COLL' ; if (PRESENT(mode_paral)) my_mode  =mode_paral

 if (Mels%nspinor == 2) MSG_WARNING("nspinor=2 not coded")

 if (PRESENT(names_list)) then
   my_nkeys=SIZE(names_list)
   ABI_MALLOC(my_keys,(my_nkeys))
   my_keys = names_list
 else
   my_nkeys=NNAMES
   ABI_MALLOC(my_keys,(NNAMES))
   my_keys = ANAMES
 end if

 ABI_DT_MALLOC(data_p,(my_nkeys))
 ABI_MALLOC(tab,(my_nkeys))
 tab=0

 my_nkeys=0; str = "  ib"; ii=4
 do ikey=1,SIZE(my_keys)
   key = my_keys(ikey)
   call my_select_melements(Mels,key,flag_p,data_p(ikey)%arr_p)
   if (flag_p==2) then
     my_nkeys=my_nkeys+1
     tab(my_nkeys)=ikey
     str(ii+1:)=" "//TRIM(tolower(key))
     ii = ii+MAX(1+LEN_TRIM(key),10)
     ABI_CHECK(ii<490,"I'm gonna SIGFAULT!")
   end if
 end do

 write(msg,'(2a)')ch10,' === Matrix Elements stored in Mels% [eV] === '
 if (PRESENT(header)) write(msg,'(4a)')ch10,' === '//TRIM(ADJUSTL(header))//' [eV] === '
 call wrtout(my_unt,msg,my_mode)
 if (Mels%nspinor == 2) then
   call wrtout(my_unt, "Sum_ab M_ab, M_11, M_22, Re(M_12), IM(Re_12)" ,my_mode)
 end if

 if (my_nkeys==0) GOTO 10
 write(fmt,'(a,i4,a)')'(1x,i3,',my_nkeys,'(1x,f9.5))' ! width of 10 chars

 do isppol=1,Mels%nsppol
   do ikibz=1,Mels%nkibz
    if (Mels%iscalc(ikibz,isppol)/=1) CYCLE

    write(msg,'(a,3es16.8,a,i2,a)')" kpt= (",Mels%kibz(:,ikibz),") spin=",isppol,":"
    call wrtout(my_unt,msg,my_mode)

    b1 = Mels%bands_idx(1,ikibz,isppol)
    b2 = Mels%bands_idx(2,ikibz,isppol)

    if (Mels%flags%only_diago==1.or.my_prtvol==0) then
      ! Print only the diagonal.
      write(msg,'(a)')str
      call wrtout(my_unt,msg,my_mode)
      do ib=b1,b2
        if (Mels%nspinor == 1) then
          write(msg,fmt)ib,(REAL(data_p(tab(ikey))%arr_p(ib,ib,ikibz,1))*Ha_eV, ikey=1,my_nkeys)
          call wrtout(my_unt,msg,my_mode)
        else
          ! Write sum_ab, then diagonal elements, finally Re_12, Im_12
          write(msg,fmt)ib,(real(sum(data_p(tab(ikey))%arr_p(ib,ib,ikibz,:)))*Ha_eV, ikey=1,my_nkeys)
          call wrtout(my_unt,msg,my_mode)
          if (my_prtvol > 0) then
            write(msg,fmt)ib,(real(data_p(tab(ikey))%arr_p(ib,ib,ikibz,1))*Ha_eV, ikey=1,my_nkeys)
            call wrtout(my_unt,msg,my_mode)
            write(msg,fmt)ib,(real(data_p(tab(ikey))%arr_p(ib,ib,ikibz,2))*Ha_eV, ikey=1,my_nkeys)
            call wrtout(my_unt,msg,my_mode)
            write(msg,fmt)ib,(real(data_p(tab(ikey))%arr_p(ib,ib,ikibz,3))*Ha_eV, ikey=1,my_nkeys)
            call wrtout(my_unt,msg,my_mode)
            write(msg,fmt)ib,(aimag(data_p(tab(ikey))%arr_p(ib,ib,ikibz,3))*Ha_eV, ikey=1,my_nkeys)
            call wrtout(my_unt,msg,my_mode)
          end if
        end if
      end do

    else
      ! Print full matrix.
      max_r = b2-b1+1
      max_c = MIN(b2-b1+1, 9)
      ABI_MALLOC(mat,(b1:b2,b1:b2))
      do ikey=1,my_nkeys
        write(msg,'(3a)')" **** Off-diagonal elements of ",TRIM(my_keys(tab(ikey)))," **** "
        call wrtout(my_unt,msg,my_mode)
        do iab=1,Mels%nspinor**2
          mat = data_p(tab(ikey))%arr_p(b1:b2,b1:b2,ikibz,iab)*Ha_eV
          call print_arr(mat,max_r,max_c,my_unt,my_mode)
        end do
        write(msg,'(a)')ch10
        call wrtout(my_unt,msg,my_mode)
      end do
      ABI_FREE(mat)
    end if

   end do !ikibz
 end do ! isppol

10 continue
 ABI_FREE(my_keys)
 ABI_DT_FREE(data_p)
 ABI_FREE(tab)

end subroutine melements_print
!!***

!----------------------------------------------------------------------

!!****f* m_melemts/melements_zero
!! NAME
!! melements_zero
!!
!! FUNCTION
!!  Set matrix elements connecting states with different irreducible representation to zero.
!!
!! INPUTS
!!   irrep_tab=Array used to select the entries that have to be set to zero.
!!     irrep_tab(ib,ik,is)=gives the index of the irreducible representation associated to state (ib,ik,is).
!!  [aname]=The name of the array to be symmetrized, by default
!!    all calculated arrays are filled.
!!
!! SIDE EFFECTS
!!  Mels= All arrays elements connecting states belonging to different irreps are set to zero.
!!
!! PARENTS
!!      sigma
!!
!! CHILDREN
!!      my_select_melements
!!
!! SOURCE

subroutine melements_zero(Mels,irrep_tab,aname)

!Arguments ------------------------------------
!scalars
 character(len=*),optional,intent(in) :: aname
 type(melements_t),intent(inout) :: Mels
!arrays
 integer,intent(in) :: irrep_tab(Mels%bmin:Mels%bmax,Mels%nkibz,Mels%nsppol)

!Local variables-------------------------------
 integer :: is,ik,ib,jb,iname,irrep_j,irrep_i
 integer,pointer :: flag_p
 character(len=NAMELEN) :: key
!arrays
 complex(dpc),ABI_CONTIGUOUS pointer :: arr_p(:,:,:,:)

! *************************************************************************

 do iname=1,NNAMES
   key = ANAMES(iname)
   if (PRESENT(aname)) then
     if (key /= aname) CYCLE
   end if

   call my_select_melements(Mels,key,flag_p,arr_p)

   if (flag_p>0) then
     do is=1,Mels%nsppol
       do ik=1,Mels%nkibz

         do jb=Mels%bmin,Mels%bmax
           irrep_j = irrep_tab(jb,ik,is)
           do ib=Mels%bmin,Mels%bmax
             irrep_i = irrep_tab(ib,ik,is)
             !
             ! Set this matrix element to zero if the irreps are known and they differ.
             if (irrep_i/=irrep_j .and. ALL((/irrep_i,irrep_j/) /=0) ) then
               !write(std_out,*)"setting to zero ",ib,jb,ik,is
               if (Mels%nspinor==2) then
                 arr_p(ib,jb,ik,is)=czero
               else
                 arr_p(ib,jb,ik,:)=czero
               end if
             end if

           end do !ib
         end do !jb

       end do !is
     end do !ik
   end if

 end do !inames

end subroutine melements_zero
!!***

!----------------------------------------------------------------------

!!****f* m_sigma/mels_get_exene_core
!! NAME
!!  mels_get_exene_core
!!
!! FUNCTION
!!  Compute exchange energy.
!!
!! INPUTS
!!  mels<melements_t>=Matrix elements.
!!  kmesh<kmesh_t>=BZ sampling.
!!  bands<band_t>=Bands with occupation factors
!!
!! PARENTS
!!
!! SOURCE

!pure function mels_get_exene_core(mels,kmesh,bands) result(ex_energy)
!
!!Arguments ------------------------------------
!!scalars
! real(dp) :: ex_energy
! type(melements_t),intent(in) :: mels
! type(kmesh_t),intent(in) :: kmesh
! type(ebands_t),intent(in) :: bands
!
!!Local variables-------------------------------
!!scalars
! integer :: ik,ib,spin
! real(dp) :: wtk,occ_bks
!
!! *************************************************************************
!
! ex_energy = zero
!
! do spin=1,mels%nsppol
!   do ik=1,mels%nkibz
!     wtk = kmesh%wt(ik)
!     do ib=mels%bmin,mels%bmax
!       occ_bks = bands%occ(ib,ik,spin)
!       if (mels%nspinor==1) then
!         ex_energy = ex_energy + half * occ_bks * wtk * mels%sxcore(ib,ib,ik,spin)
!       else
!         ex_energy = ex_energy + half * occ_bks wtk *SUM(mels%sxcore(ib,ib,ik,:))
!       end if
!     end do
!   end do
! end do
!
!end function mels_get_exene_core
!!!***

END MODULE m_melemts
!!***
